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Abstract. - We calculate analytically the phase boundary for a nonequilibrium phase tran- 
sition in a one-dimensional array of coupled, overdamped parametric harmonic oscillators in 
the limit of strong and weak spatial coupling. Our results show that the transition is reentrant 
with respect to the spatial coupling in agreement with the prediction of the mean field theory. 



Since the landmark paper of Onsager on the solution of the 2-dimensional Ising model Q , 
the search for exactly solvable models displaying equilibrium phase transitions has become one 
of the hot pursuits in statistical physics. Over the last two decades, the interest has shifted 
from equilibrium to nonequilibrium phase transitions Q. In many respects, the latter transi- 
tions are richer and more difficult to analyse. In particular, the analogue of the equilibrium 
distribution is typically not known. On the other hand, the arguments ruling out phase tran- 
sitions in one dimension for a large class of equilibrium models no longer apply, and a number 
of exact solutions for nonequilibrium phase transitions in one dimension have been found j3|. 
Our focus here is on a nonequilibrium phase transition in a one-dimensional array of coupled 
overdamped parametric oscillators. This new transition was discovered Q in the line of work 
on phase transitions induced by multiplicative noise H§|, see also [0. The phase transition 
occurs in all dimensions. Surprisingly, it was found that a mean field analysis predicts quite 
accurately the location of the phase boundary as obtained from numerical simulations for the 
one-dimensional array. Furthermore, the transition was found to be reentrant with respect 
to the coupling strength. In this letter, we derive the analytical form of the phase boundary 
in the limit of weak and strong spatial coupling. Our exact results allow for a detailed com- 
parison with the the mean field result. In particular the reentrance of the transition is firmly 
established. 

The system considered in [Q is a linear chain of overdamped particles coupled harmonically 
and parametrically excited : 

±i = [-1 + Acos(u>t + (j>i)]xi - K(2xi - x i+ i - Xi-x). (1) 

The disorder stems from a random choice of the phases <j>i in [0, 2ir}. This system exhibits 
a nonequilibrium phase transition from a quiescent phase to an exploding phase at a critical 
value of the perturbation amplitude A. In the exploding phase the amplitude of all oscillators 
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diverges with time. This is surprising because if one suppresses either the disorder (i.e., one 
replaces the (fas by a unique <fa) or the coupling (i.e., setting K = 0), no destabilization occurs. 
Hence, the instability is entirely caused by the conjugate effects of coupling and disorder. 

To analyse the model further, it is convenient to introduce a dimensionless time unit r = ujt 
and scaled variables a — A/oj and k — K/oj. Turning first to the weak coupling limit k«1, 
we switch to new variables yi, x% — y% x exp[— (1/lu + 2k)t + asin(r + $%)], and note that the 
overall time evolution of these variables: 

yi = K ( e Qsin ( T +^ + 1 )-"sm(r+0 i ) 2/i+i + e a S in(r+^-i)-a S in(T+0i) y ._ 1 ) (2) 

is slow due to the assumption that the prefactor k in the r.h.s. of @ is small. On this 
slow time scale the periodic parametric perturbations are fast and one can replace in (||) the 
following quantities by their average over one period of the oscillations: 

1 r 27T 

e Qsin(r+0 i + i)-asm(r+0 i ) __ ) x / ^ & oi sin(T+0 i+ i)-a sm(r+0;) 

2vr 







= JL f dr e 2asin ( "' +1 2 "' ) CQST = A,. (3) 



This procedure is tantamount to the Bogoliubov-Mitropolski technique, familiar from the 
literature on oscillators || . With the replacement given in (^) , (0) reduces to : 

yi = K,(\iy i+1 + Xi-iyi-i). (4) 

Clearly, for almost all initial conditions, the asymptotic time dependence of yi will be dom- 
inated by the largest eigenvalue A + of the symmetric tri-diagonal random matrix = 
Aj +1 i = Ai, Vi, the remaining elements Ajj being zero. The phase boundary for the transi- 
tion to explosion is thus given by (l/u> + 2k) /k = A + . 

In the limit of an infinite system, the value of A + is expected to be self-averaging. To 
determine this value, we note that the A^'s are statistically independent, hence our problem 
belongs to the traditional theory of random tri-diagonal symmetric matrices with independent 
random elements. Since in our case the diagonal elements of the matrix A are zero - unlike 
the famous cases of matrices coming from disordered models of linear chains with randomly 
chosen spring constants and/or masses - our matrix belongs to the "type I" disorder in Dyson's 
terminology |9|, Furthermore the support of \ is bounded, A, G [A m j n , \ m ax], and a 
theorem due to Hori and Matsuda 10, ll[ states that A + is bounded from above by 2A r 



The latter bound is realized for the choice A^ = X max ,Wi. We are however concerned with 
the typical value of A+, as it is observed with probability one in the thermodynamic limit. 
A numerical evaluation of the spectrum of a large A matrix in fact points to a value of A+ 
significantly lower than the "naive" estimation given by the upperbound (cf. fig. |l|). A 
detailed analytic analysis (cf. appendix for the major steps of the proof) however reveals 
that 2A mQX is in fact the typical value of A+ and the observed deviations in the numerics are 
entirely the result of strong finite-size effects. 

We conclude that in the limit of small k the phase boundary is given by: 

1/UJ + 2K = 2X rnax = - P dr e 2acosT = 2 Jo (2a). (5) 



i) 



7T 

This result has to be compared with the mean-field prediction |4]]: 

dre- (1/w+K)T / (2asin(r/2)) = -, (6) 
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Fig. 1 - Positive part of five realizations of 600 x 600 symmetric random matrices ; these matrices 
are of the A-type, that is tri-diagonal symmetric with zero elements on the diagonal, and random 
independent samplings in the interval [1, 2] on the off-diagonal line. The upper bound of the main 
spectrum appears to be strictly smaller than 4, in disagreement with the more careful analysis given 
in the appendix. 



which, in the region of k small and a large, reduces to: 

Io(2a) = sinh[7r(l + 1/uj)] 
y/a \/ttk 

The two asymptotic formulae are not identical, but the leading term of the asymptotic branch 
is the same, namely a ~ — ^logK. This partly explains the observed similarity of the two 
instability curves in this region. We close the discussion of the small k limit by noting that 
the above method can unfortunately not be generalized to higher dimensions : the arguments 
are very specific to dimension 1, while the variables Xi are no longer independent in higher 
dimensions. 

Next we turn to the description of the large K regime. In this case it is convenient to 
transform away the — Xi term by setting x% = yi exp(— t/uj), and to introduce the discrete 
Fourier transforms: 

i 

y,e^ N , q = 0, . . . , N - 1. (8) 

V 3=0 

The evolution equation for these Fourier modes reads : 

N — l 

Vi = a Yl yq'C q - q >(T) - (4k sin 2 y q , (9) 

q'=0 
N-l 

with C„(t) = — Y e 2lTr ^/ Jv cos(r + ^). 

In the limit of k large, the system becomes very "stiff" and any initial short wavelength 
perturbation is expected to be damped out very rapidly, cf. second term in the r.h.s. of (Q). 
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We therefore assume that, as can be verified a posteriori, the long-time behaviour of the q =/= 
modes is completely determined by their coupling to the slow zero-th mode q = 0, i.e. by 
the "reinjection" of energy through the q' — contribution in the first term of the r.h.s. of 
(^). Dropping furthermore the irrelevant contributions due to initial conditions, one concludes 
that for q 0: 

y q = a[ due-^ sia2 ^ N ^'^y (u)C q (u). (10) 



Inserting this result in the equation for q = leads to: 

N-i ,t 



= a 2 V/ due-^ s ^^/ N ^y (u)C g (u)C. q (t). (11) 

„-n Jo 



yo 

9=0 



One expects that the quantity Y^J Ne-^ Ksin ^l^-^C^C-^t) is self-averaging in 
thermodynamic limit |]l3|| , hence one can perform an average over the phase disorder directly 
in the dynamical equation: 

2 ft ( i N \ 

i0 = TJ du \jjY,^ [iKS ' n2 * q/mt ~ U) )co S (t-u)y Q (u) (12) 



a 



9=1 

2 ft / /•! 



y J du I dqe -[^m]{t-u)\ CQs(t _ u)yQ(u) _ (13) 

where the thermodynamic limit has been taken in the second equality. The resulting equation 
can be readily solved for the Laplace transform j/o^) of yo (assuming yo(t — 0) = 1) : 

Vo{z) = ^7T> (1 4 ) 

where 'g(z), the Laplace transform of g{t) = cost x f Q dq e ~[ 4Ksm2 7r '] i J is found to be 



( y/(z* + l)((z + AnY + 1) + z(z + 4«) - 1 V 
3( * )= ^ 2(z' + l)((z + 4^ + 1) ) ' (15) 

The long-time asymptotic behaviour of yo(t) turns out to be governed by the single real 
positive pole 7 of yo(z). The phase transition to explosion takes place when 7 = l/to, or more 
explicitly when 



1 a 2 a 2 Vlu-z + I + lu- 1 ]' 

- = ^-ff(V w - —7= zr— -; , (16) 

where the large n limit is invoked in the second equality. The dependence a oc k 1 / 4 explains 
the flat shape of the reentrance curve. In comparison the mean-field model predicts a be- 
haviour a ~ y2/t/w, consistent with the fact that the reentrance is expected to occur at a 
smaller value of k (for a given value of a) in the mean-field system, since the connectivity of 
the system enhances its effective stiffness. 

We finally turn to a comparison of the above derived asymptotic results with the mean 
field and numerical results, cf. fig. |[ The differences between the asymptotic ID results 
and the mean-field value are rather small. More revealing is that the ID system is less 
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Fig. 2 - Comparison between stability curves for the mean-field model (solid line), 1-d prediction for 
small (dashed line) and high (dash-dotted line) k values. The numerical results of |^| are also recalled 
(circles) (for all these curves, uj = 10). 



stable than the mean-field model for small and (very) large n values (cf. the respective 
asymptotics k 1 / 4 <C ft 1 / 2 ), while the opposite seems to be the case for intermediate values : 
if the asymptotic evaluation for k large is reliable for k ~ 4, we can see from fig. || that there 
exists a range of coupling values for which the ID system is more stable than the mean-field 
model. This inversion seems to be confirmed by the numerical simulations of the ID model 
in Q| and highlights the subtleties of the destabilisation process taking place in the system. 

We conclude that the analytic theory presented here describes the salient features of the 
phase boundary for the nonequilibrium phase transition occurring in a one-dimensional array 
of coupled, overdamped parametric harmonic oscillators with quenched random phases. The 
fact that the deviations from the mean field results are small raises the question of the critical 
dimension. Interestingly, the above presented procedure to calculate the phase boundary for 
large k can also be applied to higher dimensions. One finds a mean field behaviour for a spatial 
dimension larger or equal to 3, with the reentrance curve given by a oc a/k. The dimension 2 
is "critical" , characterized by logarithmic corrections a oc y/n/logn. This strongly suggests 
that the critical dimension of the phase transition is in fact equal to two. 

Appendix. - In this appendix, we review some spectral properties of the random matrix 
A. First, one notices that its spectrum is always symmetric with respect to zero, for if {y n } is 
an eigenvector associated with the eigenvalue A, then {(— l) n y n } is another eigenvector with 
the eigenvalue —A. We can thus restrict our attention to the negative eigenvalues —A, A > 0. 
The eigenvector equation reads — Xy n = X n Vn+i + A„_iy„_i and can be recast into the form 
z n+ i = -fj, n /(\ + z n ), with z n = A„_iy n _i/?/„, [i n = A 2 . One can now reproduce step by step 
the demonstration performed by Schmidt [l2] | , whose conclusions are that the proportion of 
negative eigenvalues comprised in [—A, 0] is equal in the thermodynamic limit to the proportion 
of negative z n in the sequence zq,z\, . . .. Turning first to the case where the /i n are constant 
(=/i) two different situations arise. For A 2 < Afi, the mapping from z n to z n+ i has no (real) 
fixed points, and as a result the positive z n are substantially visited (i.e. the proportion of 
positive z n is non zero when N — > oo). For A 2 > 4/i the appearance of a globally attractive 
negative fixed point z + = {— A + \J \ 2 — 4/i,)/2 implies that the sequence will ultimately end 
up in the R~ region. Hence, —A = is the lower bound of the spectrum, since beyond 
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this value, the proportion of negative z„-values is 1. 

We next turn to the more complicated case of A„ being independent random variables 
sampled from a distribution with support [X m i„,X max ] (with X rn i n ^ 0). Since the G 
\pmin ! MraaJ are a l so independent random variables, the sequence z n becomes a Markov pro- 
cess. In continuation of the arguments given above, the following conclusions can be drawn. 
For A > 2X max it follows that A 2 > 4[i n is true for all allowed values of fi n , and the z n are 
attracted to a region [z + (fj, max ), z + (fi min )] which lies entirely on the negative real axis. One 
concludes that the spectrum, assumed to be self- averaging in the limit N — > oo, cannot ex- 
tend beyond the value 2X max (a result also demonstrated in j|o|). Considering the situation 
where A is slightly less than 2X max , the z n will still more often than not be attracted to the 
region [— A/2, z+(/x m i n )]. Exceptionally however, when fj, n is sampled outside [/i m j„, A 2 /4], 
the z n will make an excursion toward the region ] — oo, — A]. If a long enough such series 
of successive \i n is sampled, z n will eventually cross the —A value and next become positive. 
Even though unlikely, the probability for these events is not zero in thermodynamic limit, for 
their occurrence only requires a finite number of steps in the sequence of z n . As a result, one 
concludes that the effective boundary for the spectrum is indeed ±2X max . Moreover, one can 
understand qualitatively why the spectra of figure ^ do not seem to saturate this bound as 
follows. The number of successive steps with /j, n G I = [A 2 /4, /j, max ], required to reach R + 
can be considered more or less as constant (with respect to A), say Uq. Then the probability 
for observing a positive value of z n is of order p n ° , where p is the probability to sample \x n in 
/. In the example of figure |l|, p is simply given by (X max - A/2)/(A 

max ^min 

) = 2 - A/2. 

Moreover, a rough evaluation of no gives at least 4 steps. Then the proportion of eigenmodes 
greater than 3.5 can be (very generously !) estimated to be (2 — 3.5/2) 4 = 0.3% . . ., explaining 
the apparent absence of these eigenvalues in figure [jL 
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